Quantum Simulation of the Hubbard Model: The Attractive Route 
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We study the conditions under which, using a canonical transformation, the phases sought after for 
the repulsive Hubbard model, namely a Mott insulator in the paramagnetic and anti-ferromagnetic 
phases, and a putative d-wave superfluid can be deduced from observations in an optical lattice 
loaded with a spin-imbalanced ultra-cold Fermi gas with attractive interactions, thus realizing the 
attractive Hubbard model. We argue that the Mott insulator and antiferromagnetic phase of the 
repulsive Hubbard model are easier to observe in the attractive Hubbard mode as a band insulator 
of Cooper pairs and superfluid phase, respectively. The putative d-wave superfluid phase of the 
repulsive Hubbard model doped away from half-flUing is related to a d-wave antiferromagnetic 
phase for the attractive Hubbard model. We discuss the advantages of this approach to 'quantum 
simulate' the Hubbard model in an optical lattice over the simulation of the doped Hubbard model 
in the repulsive regime. We also point out a number of technical difficulties of the proposed approach 
and, in some cases, suggest possible solutions. 



PACS numbers: 

I. INTRODUCTION 

Understanding the phase diagram of the two- 
dimensional (2D) single-band Hubbard model is consid- 
ered by many as the 'Holy Grail' of the theory of strongly 
correlated systems. In the most interesting regime, this 
model describes a system of spin-^ {i.e. two species of) 
fermions hopping on a 2D square lattice with repulsive 
(on-site) interactions, and average lattice filling less than 
one fermion per site. This model has been proposed as 
the minimal model that explains the observation of d- 
wave superconductivity with fairly high critical temper- 
ature in the doped cuprate materials [1] (for a review on 
doped Mott insulator, see 0). At half filling, where only 
one particle per site is allowed, it is by now rather well 
established that the model is a Mott insulator, which at 
low temperatures (below a characteristic scale, the Neel 
temperature, TNeei) orders anti-ferromagnetically. Away 
from half filling, the nature of the ground state is a sub- 
ject of heavy debate. One of the most challenging open 
issues is whether the Hubbard model on a 2D square lat- 
tice would support a d-wave superconducting phase at a 
(relatively) high temperature. The fact that the Hubbard 
model can support such an instability has been theoreti- 
cally proven in double-chain systems coupled by hopping 
(known as two-leg Hubbard ladders) [31 HI [51 [3] . However, 
whether this result extends to the 2D model consisting 
of an infinite number of coupled chains is still extremely 
controversial [H [71 [5] . At present , neither analytical nor 
numerical studies are able to settle the issue. 

Due to the spectacular advances in the optical ma- 
nipulation of ultra-cold atomic gases, one very promis- 
ing route for studying the low temperature phases of 
the Hubbard model has opened up recently [F]. Indeed, 



ultra-cold Fermi gases loaded into an optical lattice can 
be regarded as almost ideal quantum simulators of the 
Hubbard model, where independent control of the hop- 
ping amplitude, t, and the on-site interaction energy, U, 
are both experimentally available. Exploiting this fact, 
the Mott insulating phase of the Hubbard model has re- 
cently been demonstrated in a 3D cubic optical lattice 
(where the center is at half-filling) by several experimen- 
tal groups [inilTT]. Many other groups are currently en- 
gaged in similar experimental endeavors with the 
main focus on realizing the repulsive Hubbard model on 
a 2D square lattice away from half- filling, namely the 
regime where d-wave superfluidity is speculated to exist. 
However, one of the main problems that lie ahead in this 
program has to do with the currently accessible temper- 
atures for the Fermi gases in optical lattices. At present, 
these temperatures (of the order of a few tenths of the 
Fermi energy of a non-interacting gas of similar average 
density) still largely exceed the Neel temperature (Tnooi), 
thereby washing out any anti-ferromagnetic order in the 
half-filled system [13]. 

However, there are other problems with the present ap- 
proach to simulate the repulsive Hubbard model, which 
seem not to have received so much attention thus far. 
One of the most remarkable ones is the difficulty of dop- 
ing away from half-filling a Fermi gas loaded in an op- 
tical lattice. In this case, the situation is very different 
from doping in solids, mainly because of two reasons: 
i) The existence of an overall harmonic trapping poten- 
tial superimposed on the optical lattice potential that 
makes the system inhomogeneous and tends to favor the 
maximum site occupancy (i.e. two fermions per site) 
near the center of the trap. As recently demonstrated 
experimentally [TDl 111] , for small number of particles this 
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tendency can be balanced by the onsite repulsion energy, 
which yields a Mott insulator near the center surrounded 
by a 'metallic' region where the density of holes is non- 
uniform, ii) Although the number of available lattice 
sites can be controlled with accuracy, the total number 
of atoms is still hard to measure accurately and it is also 
subject to variations from shot to shot that are inherent 
to the preparation process [14 . 

Another source of problems has to do with the need to 
independently control the on-site interaction, U , and the 
hopping amplitude, t. As the Neel temperature (Tnooi) 
(below which the system orders anti-fcrromagnetically, 
and, upon doping, the putative d-wave superfluid may 
appear), scales as TnccI ~ t^/U, and thus rapidly de- 
creases if the ratio t/U is made very small by increasing 
the optical lattice depth, it is desirable to have indepen- 
dent control of both t and U . In order to achieve this, the 
s-wave scattering length, Os {U (x Qs, roughly speaking) 
that characterizes the strength and sign of the atom-atom 
interaction, must be tuned towards a Feshbach resonance 
where as — > ±00. Since the current interest is in real- 
izing a Hubbard model with repulsive interactions, the 
side of the Feshbach resonance where the atom- atom in- 
teraction is repulsive (i.e. > 0) must be used (see, 
however. Sect. Ill for further remarks). On this side of 
the resonance, there is a weakly bound molecular bound 
state |15j . with which the atoms in the continuum (that 
is, in the lowest and highest Bloch bands, when loaded in 
a lattice) have a sizable overlap near the resonance. Thus, 
at sufficiently low temperatures, Feshbach molecules form 
resulting from three atom collisions [15]. A collision of 
one of these molecules with a third atom can cause the 
molecule to make a transition into a more bound molec- 
ular vibrational state. The released energy is taken away 
by the colliding atom, which therefore causes undesir- 
able heating of the system. Also, the presence of these 
molecules is not accounted for by the single-band Hub- 
bard model, which is the goal of the quantum simula- 
tion. Furthermore, as the scattering length increases, U 
also increases and become of the order of the separation 
between Bloch bands, thus leading to the break-down of 
the single band approximation [16l [171 [18]. 

In this article, we propose to explore a different route 
to simulate the Hubbard model in a regime where the 
onsite interaction is attractive. Theoretically, the attrac- 
tive and repulsive regimes are related by a transformation 
that is well known in the literature of the Hubbard model 
and it is, for completeness, reviewed in Sect.[llj More re- 
cently, in the context of cold atomic gases, this transfor- 
mation has been used by Moreo and Scalapino [19j , who 
pointed out the connection between the Fulde-Ferrel- 
Larkin-Ovchinikov state in the attractive Hubbard model 
and a state with stripes in the repulsive model. These 
authors also briefly considered the relationship between 
the d-wave superfluid order parameter for the repulsive 
model and a d-wave anti-ferromagnetic order in the at- 
tractive case. For the one-dimensional Hubbard model 
(note, however, that there is no d-wave superfluid phase 



in this case) , the transformation was also used in Ref. f5D] 
in an analysis of the noise correlations of the attrac- 
tive Hubbard model with spin imbalance. Moreover, 
the physics of the attractive Hubbard model has also at- 
tracted much interest by itself j21j, especially in recent 
times and in connection with cold atomic gases and the 
physics of the BEG to BCS crossover [22] . 



In this article, we explore in depth the possibilities of- 
fered by the attractive model to understand the physics 
of the repulsive Hubbard model. We pay special atten- 
tion to the effects of the trapping potential as well as the 
peculiarities of the physical realization of the negative U 
Hubbard model in optical lattices. We thus argue that 
the attractive regime presents a number of advantages 
for the quantum simulation of the Hubbard model in an 
optical lattice. We also discuss how the (negative-C/) 
equivalent phases of the paramagnetic Mott insulator, 
the anti-ferromagnetically ordered Mott insulator, and 
the putative d-wave superfluid may be observed. 



We would like to emphasize that one of the main ad- 
vantages of the route suggested here is that the equiv- 
alence of doping away from half-filling the optical lat- 
tice system in the repulsive regime can be achieved by 
creating a spin-imbalanced gas. Fermi gases with spin- 
imbalanced populations are nowadays routinely created 
in the laboratory |23j . and the magnetization (which is 
fixed for the duration of the experiment by the prepara- 
tion method) can be controlled to within a few percent 
accuracy. Other advantages will be discussed further be- 
low. However, our approach does not solve the problem 
of how to achieve lower temperatures for the fermions on 
the lattice Uni. 



The outline of this article is as follows. In Sect. [TTj 
we discuss the transformation that formally maps the re- 
pulsive Hubbard model into an attractive one. We also 
define the two attractive Hubbard models we shall be 
concerned with in this paper, together with their corre- 
sponding repulsive models. We also describe how vari- 
ous physical operators and order parameters are affected 
by the mapping. Some important caveats concerning 
the realization of the attractive Hubbard model are dis- 
cussed in Sect. |III[ In Sect. jlVj we discuss the equiva- 
lent state of the paramagnetic phase of the Mott insu- 
lator (as well as possible ways of detecting it), whereas 
in Sect. [V] we do the same for the equivalent state of 
the anti-ferromagnetically ordered Mott insulator. The 
effect of doping with holes, which, as mentioned above, 
amounts to a spin-imbalanced situation in the attractive 
case, is analyzed in Sect. |VI| Finally, in Scct. |VIll we offer 
the conclusions of the present work as well as mentioning 
some open problems. 



3 



II. THE HUBBARD MODEL AND THE 
PARTICLE-HOLE TRANSFORMATION 

The Hamiltonian of the single-band Hubbard model 
reads: 

H = -tY. 4c,-+C/^("^T-^)(^a-^)+^^ext, (1) 

(ij) o i 

where t is the hopping amplitude and U the one-site in- 
teraction. We consider here only the case where the sites 
i of the lattice constitute an hypercubic lattice (square in 
two dimensions, and cubic in three dimensions); in 
the hopping term means that the sum runs over nearest 
neighbor sites only. In the above equation, riicr = cl^Cia 
is the occupancy of spin a =1 , | fermions at the «-th site. 
For simplicity and unless otherwise stated, we will refer 
in the following to the two dimensional case, and thus 
to a square lattice. All our results are straightforwardly 
generalizable to the case of any hypercubic lattice. 

We have denoted as i/cxt all the external fields like 
chemical and trap potentials as well as an external Zee- 
man (magnetic) field that act upon the system [41]. Their 
effects will be discussed below. In the grand canonical 
ensemble, 

Fext = ~t^N + hM (2) 

where the total number operator N — ^ and the 
total magnetization, M = and /i along 

with h are determined by the condition that averages of 
N and M over the grand canonical ensemble yield the ex- 
perimentally observed values |42j . However, cold atomic 
gases are prepared in eigenstates of both N and M, and 
therefore the relevant ensemble is the canonical instead 
of the grand canonical. Although it is important to keep 
this distinction in mind, we expect that for sufficiently 
large N, the results of both ensembles coincide and thus, 
we shall use the grand canonical for the calculation of 
the experimental signatures of the different phases to be 
described below in Sect. IVII 

The Hamiltonian of Eq. ([T|) has been written in a form 
such that for a uniform system the ground state will have 
exactly one particle per site, at any temperature, for U > 
and i?cxt — 0. However, note that in real experiments 
cold-atomic gases are harmonically trapped. Therefore, 
the most general form of i?ext reads: 

-ffext =^{ei-^i){ni^+nii-l)-^hi{ni^ - n^i) , (3) 

i i 

where e; is the shift in the local chemical potential caused 
by the trap. We have added an unimportant constant to 
the total energy ( — ~ A*) ); which will become 

convenient further below. In current optical lattice ex- 
periments, we have = \eoi'^, but more general forms 
of the trap may become available in the future. In the 
case relevant to experiments, the Zeeman field hi = h is 
uniform (and it is used to adjust the total magnetization 



in the grand canonical ensemble). However, in Eq. ^ 
we have assumed it to be site-dependent for further con- 
venience. 

We next note that, formally, the sign of the interac- 
tion term (the term cx C/ in Eq. [l]) can be changed by 
means of the following (particle-hole) transformation on 
a bipartite lattice such as the 2D square lattice: 

Ci| < > Ci^. 

Note that the transformation leaves the operators of the 
spin I fermions unchanged. However, it affects the spin | 
occupation operator: n^i <-^- 1 — n^, and thus the sign of 
the interaction term changes {U ^ —U) while hopping 
term is left invariant (the minus sign in the right-hand 
side of Eq. ^ takes care of this). However, it will be 
important for the discussion that follows that the trans- 
formation exchanges the roles of hi and {ti — ^) in Eq. (|3|. 
Mathematically, 

^ (5) 
-2_^hi (n,| + n^i - 1) . 

i 

If we insist in using the point of view of the canonical en- 
semble, the transformation ^ implies that N ^ N' = 
M + Af and M A/' ^ N - Af, where TV is the total 
number of lattice sites. Thus, if we consider an unmagne- 
tized {i.e. M = 0) system where TV = (1 -I- x)Af {x being 
the doping, where a; = corresponds to half-filling in the 
uniform case), we have that N' = Af and M' = xAf. In 
words, the doped lattice at J7 > away from half-filling 
maps onto a, U < system at finite magnetization. We 
emphasize that, as discussed below, the details of the or- 
der that the system develops depend not just on usual 
factors such as the temperature and strength of the trap- 
ping and lattice potential, but are also constrained by 
these globally conserved quantities. 

It is also convenient to recall that, in momentum space, 
the transformation of Eq. Q becomes: 

Cfci ^ 4+Qi ' "^fcT ^ Cfc| , (6) 

where Q = (tt/o, tt/u) is the nesting vector with a the lat- 
tice spacing. This expression can be used to obtain the 
way the different order parameters and the correspond- 
ing phases transform between the U > and the U < 
cases. This is shown in table |Tj The way the transforma- 
tion affects the different phases expected for the Hubbard 
model is also illustrated in Fig. [l] 

Finally, for the sake of clarity, we will spell out in what 
follows the two physically different attractive Hubbard 
models considered in this article, as well as the repulsive 
Hubbard models onto which, via the particle-hole trans- 
formation Eq. [4] they are mapped. The first attractive 
Hubbard model (from here on called model Al) has the 



following Hamiltonian (in the grand canonical ensemble) : 

(ij) a i 



-^(ej - ^l) (n^i +714 - 1) 

i 



(7) 



This is the Hamiltonian that describes cold atomic sys- 
tems trapped in an optical lattice, with the overall har- 
monic trapping potential = ^egi^, and a uniform Zee- 
man field that can be viewed as a knob to tune the spin 



imbalance (See Sect. Ill I. In this article, we argue that 
simulating this Hamiltonian has a number of advantages 
over current attempts at simulating the repulsive Hub- 
bard model in presence of the harmonic trap (from here 
on called model R2), whose Hamiltonian reads: 



(ei - ^i) (n,| -f - 1) 



(8) 



The other attractive Hubbard model (from here on called 
model A2) is described by the following Hamiltonian: 



(9) 



This is the Hamiltonian that can be obtained from the 
repulsive model R2 (Eq. [s]) via the particle-hole trans- 
formation of Eq. ^ . Model A2 has an inhomogeneous 
Zeeman field stemming (via the transformation) from the 
trapping potential of model R2. In this article, model 
A2 is used to help us understand e.g. the coexistence of 
phases in the current experimental regime of the repul- 
sive Hubbard model (R2). Note that the attractive model 
Al that we are advocating does not map onto the cur- 
rently studied repulsive Hubbard model R2, but instead 
onto the repulsive Hubbard model in a inhomogeneous 
Zeeman field (called model Rl here) : 



{ij) a i 

+ Yi<^i- M) ("•jT - "U) 
i 

-hY.{n,^ + nn-l). (10) 

i 

Thus, to summarize, using the mathematical transfor- 
mation of Eq. Q , we can relate the physics of the attrac- 
tive Hubbard models to that of the repulsive ones. In 
particular the observation of one particular phase {e.g. 



(a) 



(b) 



REPULSIVE HUBBARD MODEL ATTRACTIVE HUBBARD MODEL 

ParsnagiuUc Motl Insiilatgr ( T > ) Disordmd CDW ( T > T,"^™' ) 

Anti-fcrmniagiiclic MoH Insulator ( T < T^^^, ) Orditred CDW ( T < T^*^"'^' > 

cc> §\ /\8/ 

Dopf^d Ajitj-fcrraniagiictic Mott InsuJatoi' 
(d-wavc supcrfiuLd) 



Spin polarized CDW{d-wavc aiiti-fcrroiiiagnct) 



FIG. 1: (Color online) Schematic diagram of how the trans- 
formation of Eq. Q works. In (a) we sketch how it acts on 
a single site {e.g. red stands for spin up and blue for spin 
down), and in (b), the way it transforms different types of 
states on a uniform optical lattice. In addition to the phases 
depicted above, an anti-ferromagnetic state at f/ > ordered 
along the x or y direction corresponds to a superfluid phase 
of fermion pairs. The bottom diagram illustrates one of the 
main points of this paper, namely that doping the attractive 
Hubbard model corresponds to introducing spin inbalance in 
the repulsive Hubbard model. 



a d-wave AF phase) at 17 < would directly imply the 
existence of the corresponding phase at C/ > (the pu- 
tative c?-wave superfluid phase). Below we shall see that 
the realization of some of these phases in the attractive 
regime requires sometimes less stringent conditions than 
the corresponding ones in the repulsive regime. Further- 
more, as described above, the exchange of roles of the 
Zeeman field and the chemical potential terms effected 
by the transformation implies that we can simulate the 
doping of the Mott insulator by creating a system with 
a finite magnetization {i.e. a spin-imbalanced system). 
Also, even if the temperatures that can be achieved in 
current experiments do not allow for the investigation 
of the low-temperature ordered states, we may expect 
that, by the attractive route, some useful insights can 
be gained into other controversial issues for the high- Tc 
community, such as the nature of the normal state of the 
2D repulsive Hubbard model away from half-filling. 



III. REALIZATION OF THE ATTRACTIVE 
HUBBARD MODEL 

In principle, since the onsite interaction energy U is 
naively proportional to the atomic scattering length in 
free space [TS], Us, the U < regime can be accessed by 
sweeping the magnetic field to the side of an inter-species 
Feshbach resonance jTS] where Os < 0. In the literature 
of the BCS to EEC crossover [fsi [H il], this side is 
known as the "BCS side" of the resonance. 
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TABLE I; Transformation between phases (and their asso- 
ciated order parameters) of the Hubbard model for U > 
and U < under the canonical transformation of Eq. Q 
(cf. Eq. [g] in momentum space) . In the expressions below, 

(l>lf^ = (cOSfco; 



COS ky) is the d^2_y2 lattice form factor. 



Feshbach molecules [TS] leads to inelastic losses. 

To summarize, the attractive regime can be reached by 
making the scattering length as negative, but not beyond 
a certain limit where for Qg = —I*, {I* depending on the 
lattice dimensionality and other parameters p5 j [26 | l27]) 
the scattering amplitude has a geometric resonance. 



u >0 


u <0 


paramagnetic Mott 
insulator (MI, rii = 1) 


disordered lattice with 
rii — or rii — 2 


s-wave AF MI (z) 


s-CDW insulator 

A (s)CDW / t \ 


s-wave AF MI (x,y) 


S-wave superfluid 


d-wave superfluid (?) 


d-wave AF metal (?) 



However, the above point of view entirely neglects the 
subtleties of the scattering problem on a lattice poten- 
tial, as it turns out that U is not a linear function of 
tts, in the general case [251 US [23 • The details of the 
dependence of the zero momentum scattering amplitude, 
fo{as) (and U cx /g), on the atomic scattering length, 
as, are determined by the dimensionality and other pa- 
rameters of the lattice [251 (23 [2Z]- However, all these 
results share one common feature, namely the existence 
of a particular length scale, > 0, such that for — — Z* 
the scattering amplitude /o exhibits a (geometric) reso- 
nance [251|2ni[lZ]: /o oc as/{l + Qs/l*)- Indeed, the res- 
onance can be approached either by changing the lattice 
parameters, or by changing through a Feshbach reso- 
nance. We shall focus on the latter case here. To realize 
the attractive Hubbard model, U < 0, we require tun- 
ing the scattering length to the attractive side (os < 0) 
but such that \as\ < h. As we approach the Feshbach 
resonance and /o diverges to — cxd, crossing the geometric 
resonance beyond leads to the interaction becoming 
effectively repulsive, and as a consequence, close to the 
geometric resonance, a weakly bound state appears. The 
existence of this bound state has been discussed in the 
literature [251 ES], and if the temperature is sufficiently 
low (compared to the binding energy of the bound state), 
many bound states will be created even if the effective 
interaction between the atoms in the lattice is repulsive. 
This regime is clearly not described by the repulsive Hub- 
bard model, because it does not take into account the 
bound states. The situation is similar (although probably 
less harmful for the system [25 ) to the one encountered 
as Qs ^ +00, on the so-called BEC side of the Fesh- 
bach resonance. In this regime, the scattering amplitude 
corresponds to that of a repulsive effective interaction, 
which may lead us to think that the system is described 
by the repulsive Hubbard model, except crucially, for the 
existence of the lattice molecular bound states. But in- 
deed, as described in the introduction, the existence of 



IV. ANALOG OF THE MOTT INSULATOR IN 
THE ATTRACTIVE REGIME 

Let us first start by looking at the U < system with 
a balanced population of spin up and down fermions 
(Model AI, Eq. [7]with h = 0) and at temperatures T 
where \U\ can be made such that T < \U\ (however, 
in this section T is assumed to be large compared with 
t'^/\U\). Using the transformation, Eq. (|4|, this corre- 
sponds to a half- filled system for J7 > (model Rl, Eq. 10 
with h — 0). This is the situation where it is known that 
a Mott insulator appears in model Rl. For U > the ex- 
istence of the Mott insulator means that the states with 
more than one particle per site are strongly disfavored 
due to the large on-site repulsion, U. The correspond- 
ing situation for /7 < is that the only allowed states 
for every lattice site are either zero or doubly occupied 
states, as shown in Fig|l] The existence of the Mott in- 
sulator in the repulsive regime thus corresponds in the 
attractive regime to having all fermions form pairs [43]. 
For sufficiently large \U\ these pairs are tightly bound, 
which means that their existence can be probed by send- 
ing photons to photo-associate them into dimers, which 
are no longer trapped, and therefore can be detected as 
a loss of atoms from the lattice [25]. 

It is worth noticing that this system of pairs exhibits 
a pairing gap (~ \U\, for large \U\) to spin excitations. 
Therefore, a measurement of the single particle proper- 
tics, like the single-particle spectral function which is 
accessible by photoemission-like spectroscopy proposed 
in [29 and recently applied to ultra-cold Fermi gases by 
the JILA group [30 , should be able to detect the pairing 
gap. However, since in model Al the harmonic trap only 
couples to the atomic density (cf. Eq. [t]), it does not 
lead to the breaking of the pairs and, therefore, it does 
not take the system out of the subspace where rii = 
or Hi — 2 (which, by virtue of Eq. |4] corresponds to half- 
filling, i.e. Tii = 1). Indeed, the role of the trap is to lift 
the large energy degeneracy (for t = 0) oi the states in 
this subspace. In other words, in absence of the trapping 
potential in model Al, and since the chemical potential 
fi couples to the total particle number, all states with the 
same total number of fermions and only doubly occupied 
or empty sites are degenerate. The trap breaks this de- 
generacy and selects as the ground state of model Al the 
state where all pairs uniformly occupy the lattice sites at 
the trap bottom. In the strong coupling limit, this state 
can be regarded as a 'band insulator' of the pairs (see 
Fig.|3|. 

A complementary way of arriving at the same conclu- 
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sion relies on the transformation of Eq. ([4|. After the 
transformation, model Al at /i = is mapped to model 
Rl , that is a repulsive Hubbard model but in a inhomo- 
geneous Zccman field (note that h ~ Q in this model too, 
but in this case it couples to the total density) . Although 
the Zeeman field affects the magnetic ordering by ferro- 
magnetically polarizing the fermions at the center of the 
trap, it does not lead to doubly occupied sites, and there- 
fore it does not affect the characteristic incompressibility 
of the Mott insulator, which depends on the existence of 
a gap [/, for large U) to all density excitations. 

It is worth to contrast the situation described in pre- 
vious paragraphs with the one found in current exper- 
iments, which are performed in the repulsive regime of 
the Hubbard model (model R2, Eq. |8|. In such a sys- 
tem, one needs to adjust the chemical potential /i (that 
is, the number of fermions, N) in order to have a half- 
filled lattice with one fermion per site at the trap bottom. 
Otherwise, at too large a n/U {i.e. large iV), the system 
energetically prefers to pay the energy cost of having dou- 
bly occupied sites near the bottom of the trap, rather 
than accumulating them far from the center, where the 
trapping energy is very large. Thus, the lattice at the 
center of the trap ceases to be in the Mott insulating 
phase, becoming a band insulator. On the other hand, 
for /i too small (i.e. small N) the optical lattice is not 
uniformly occupied by one fermion site. Testing for the 
existence of the Mott insulator at L/ > thus requires 
checking for the absence of doubly occupied sites, which 
has been already achieved experimentally by the Zurich 
and Mainz groups [101 E]. However, testing the absence 
of holes (which may appear due to thermal or quantum 
fluctuations, especially as iV or f/ decrase) is a more 
difficult task. In this regard, in the attractive regime, 
only the absence of singly occupied states needs to be 
tested. Such measurement, which implies recording the 
spatial distribution of lattice sites with different occupa- 
tions may be accessible through spectroscopic techniques 
similar to those employed to observe the 'wedding-cake' 
structure of the Bosonic Mott insulator [311 132] • 

In Figs. 3J2 we show illustrative plots (in the large \ U\/t 
limit and for T < \U\) of the experimentally measur- 
able/measured density profiles n(r) = n|(r)-|-n|(r), both 
in the attractive regime with no spin imbalance (Fig. [s] 
for model Al at h = 0) and in the repulsive regime, the 
half-filled lattice (Fig. [2] for model R2). [44 We have 
used data similar to those in current experiments, e.g. 
the Zurich group [10] (see figure captions). 

It is noticeable that for U > (model R2) , for the given 
number of particles vs. trap energy cq = muP'c?' (m is 
the mass of a single fermion and ld is the trapping fre- 
quency), there is single occupancy in the central region, 
which is the Mott insulator region. At temperatures that 
are low compared to [/, but higher than any second or- 
der (exchange, etc.) processes at the energy scale ~ jU 
(i.e. at currently achievable temperatures), the t and J, 
fermions are equally likely to be found on any site of the 
Mott-insulating region, as there is no magnetic ordering. 




80 a 



FIG. 2: Density profile for the repulsive Hubbard model in 
a harmonic trap (model R2), as a function of the radial dis- 
tance from the trap center r in units of a, the optical lattice 
spacing. Dashed line: T/[/=0.05, sofid line: T/U = 0.1, dot- 
dashed line: T/U = 0.2. Total number of particles is N=6599. 
Trapping energy eo/2U = 0.0003. 



The latter can only emerge as the temperature is lowered 
below the Neel scale, T^eei ~ t^/U. As we move towards 
the edge of the sample, the site occupancy deviates no- 
ticeably from one fermion per site, and a "metallic" shell 
appears. Its width depends on the temperature, as well 
as t,U and the trapping energy. 

On the other hand, for the [/ < case of model Al at 
/i = 0, the center of the lattice is filled with fermion pairs 
resulting from the attractive on-site interaction. Thus, 
in model Al (Fig. [S]) the effect of temperature (in the 
limit where one can disregard the second order effects 
due to the hopping) , is to create a finite density of empty 
sites near the edge of the 'band insulator' of pairs. This 
explains the deviation of the site occupancy from nir) — 
2 observed in Fig. |3] as the distance to the center, r, 
increases. We also notice that, because the trap is not 
in competition with the pairing gap as described above, 
the stability of the band insulator state is not threatened 
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FIG. 3: Density profile for the attractive Hubbard model in a 
harmonic trap (model Al at ft = 0), as a function of the radial 
distance from the trap centre. Dashed line: r/l7=0.05, solid 
line: T/U = 0.1, dot-dashed line: T/U = 0.2. Total number 
of particles is N=11728. Trapping energy to/2U = 0.0003. 
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by the trap (contrary to the Mott insulator a,t U > 0), 
and thus the fermion numbers in this regime (U < 0) can 
be larger than those at ?7 > 0, which is also illustrated 
by the sizes used to generate the figures. However, for 
model R2, the trap, as described above, tends to favor 
the band insulator in the middle for large enough N. 
When the total number of fermions N becomes smaller 
than a critical value (see next section), consideration of 
the effects of hopping will be required. 



a) ® 



spin into paper {-z) 
spin out of paper (+z) 



b) 




empty 




V. ANTIFERROMAGNETIC ORDER FOR U >0 
AND ITS ANALOG IN ATTRACTIVE CASE 

As discussed in the previous section, for large on-site 
attraction, the only two possible states for a single site 
are empty sites and doubly occupied. Singly occupied 
sites are separated from these states by the (large) energy 
gap ~ \U\. Thus all the fermions are paired and can be 
regarded as hard core bosonic entities hopping from site 
to site with amplitude ^ t^/\U\. To see this, recall that, 
in the repulsive case |33| and ignoring for the moment 
the overall harmonic trapping, the low energy effective 
model for the U > half-filled Hubbard model (with 
equal proportion of the two species) is a spin-1/2 nearest- 
neighbor Heisenberg model H^siU > 0) = JJ2{ij) ' 

Sj, J = At^/\U\. This model transforms for [/'< 0, ter A^^* j^^pg ^ 



FIG. 4: (Color online) Schematic diagram of the shell struc- 
ture of a) the half-filled repulsive Hubbard model in a site- 
dependent Zeeman field (model Rl), and b) the corresponding 
attractive model Al state configuration in a trap. 



state and the SF state for the U < case, which one can 
see formally from the order parameter mapping, using the 
transformation Eq.(6]): the checkerboard order parame- 
hc antifcrromagnetic order parameter 



upon applying Q, into a half-filled lattice described by 
the following effective Hamiltonian in terms of hard-core 
bosons {bi = Ci|Ci|): 



HMU < 0) 



r3 



,(11) 



where means that the sum runs over nearest neigh- 
bor sites only. For the U > case, the ground state 
of the Heisenberg model is a (s-wave) antiferromagnet 
(sAF). In the absence of any terms in the Hamiltonian 
that distinguish spin up and down fermions, that is, for a 
spin-isotropic Hamiltonian, the staggered magnetization 
can point in any of the spin directions, x,y, or z. 

For the U < case, the resulting system, Eq. (Ill, 



is thus a system of hard-core bosons (arising from the 
pairing of two different spin fermions) hopping on the 
lattice with the kinetic energy ~ J, and experiencing 
nearest neighbor repulsion interaction also ^ J. Several 
phases can be realized in such a system. The resulting 
nearest neighbor repulsion and kinetic energy in Eq. (11) 



in the z spin direction M^^g^g, while the SF order parame- 
ter A*^''', A'^*)^^ maps to the AF order along — , + direction, 

^stag = ^ftag + ^Afs'tag and M.^.g = M^.g - iM^,,^ m 
See Table H 

However, the presence of the harmonic trap on the at- 
tractive side leads to interesting effects. As can be seen 
from Eq. ([s]), in model Al Eq. ([T]), the trap acts as a 
chemical potential for the pairs in the U < case, and 
thus favors a completely filled center of pairs (which is 
just a band insulator in the center, cf. Fig.|3|, rather than 
either a CDW or superfluid state. In fact, in the limit 
where tunneling between sites is suppressed, the ground 
state in the trap potential is this band insulator, and only 
the kinetic energy and the nearest neighbor repulsion of 
order J prevent this state to occur. Mapping back to 
the U > case, one sees that the trap transforms into 
a Zeeman field along the z direction to become model 
Rl Eq. (10 1. This Zeeman field will lift the degeneracy 



favors an alternating pattern of empty and doubly occu- 
pied states, i.e., the checkerboard state in 2D also known 
as the "charge density wave" (CDW) state. This alter- 
nating pattern of empty and doubly occupied states for 
[/ < is the analogue of the antifcrromagnetic (along z- 
direction) state for ?7 > as shown in Figjljand Table |l] 
Another alternative ground state is simply an s-wave su- 
perfiuid of these bosonic pairs (sSF). In fact, the degener- 
acy for ordering along any direction x,y,z of the sAF for 
the [/ > case maps to a degeneracy between the CDW 



between the various magnetic states and then polarize 
ferromagnetically the spins rather than favor an antifcr- 
romagnetic order. This competition is summarized in 
Fig.i 

This effective Zeeman field in model Rl goes from large 
and positive in the center of the trap to large and nega- 
tive in the periphery. In the center of the trap one would 
thus have all the spins polarized up. This phase cor- 
responds via the transformation, to a pair on each site 
and thus to the band insulator of pairs in model Al. 
Whether the effective field in the center of the trap is 
enough to polarize fully the spin depends on the num- 
ber of particles in the trap and will be discussed below. 
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Further from the center, the Zeeman field decreases and 
becomes negative. One has thus in a certain radius a 
shell where the spins are not fully polarized. In this re- 
gion the spins preserve an antiferromagnetic order in the 
direction perpendicular to the Zeeman field albeit with a 
reduced amplitude {i.e. the AF order parameters Mq^^ 
are non-zero). In other words, in this shell, the system of 
model Rl thus possesses antiferromagnetic order in the 
X, y direction in spin space. After the transformation this 
X, y-antiferromagnetic phase for U > maps for J7 < to 
an s-wave superfiuid phase, as shown in Fig. |4] We thus 
see that looking for antiferromagnetism in the repulsive 
Hubbard model amounts to probing for superfluidity in 
the U < one. Then beyond a certain radius the trap 
prevents the pair to exist, and the system is empty (see 
Fig. [3|. In the U > system, this corresponds to a Zee- 
man field that is large enough to fully polarize to down 
the system of spins [46] . 

In summary, on the U < side the trap does not really 
spoil the search for the analogue of the antiferromagnetic 
phase. The corresponding phase is now a simple s-wave 
superfiuid that can be probed by similar techniques that 
have already been used in the continuum. The extra 
phases induced by the trap only potentially reduce the 
spatial extent of the superfiuid shell, but since they both 
corresponds to insulating regions (either band insulator 
or empty region) , they should not spoil the observation of 
the superfiuidity. To minimize the central band insulator 
region for the attractive model Al (and hence maximise 
the superfiuid signal), the chemical potential at the trap 
center must not be too strong, so that when tunneling is 
turned on, the resulting pair hopping and pair repulsion 
energy scale J (see Eq. 1 1 1 can overcome the trapping 
potential to delocalize the band insulator. First, in the 
strong coupling limit, consider the case where the whole 
trap consists of mostly the band insulator. We can es- 
timate the instability of the band insulator to this scale 
J: taking a pair at the edge of the band insulator system 
{R = R2) to a, site just beyond the edge (by a lattice 
spacing a) costs an energy ^ eo(i?2/a), and this is to be 
balanced against ,/ — 4t^/|[/|. For a large enough sys- 
tem, 7r(i?2/a)^ = N/2 where N/2 is the total number of 
pairs. Thus the band insulator s tate w ill become unsta- 
ble when roughly, t^/|C/| > e^^/N /2tt. However, for a 
full SF state in the central region of the trap, one needs 
to move outward more than just the pairs at the bound- 
ary of the band insulator. A similar estimation indicates 
that the V^iV in the above criterion is replaced by linear 
in N. Thus for a given trap energy eg, there is an up- 
per limit to total number of fermions that can be loaded 
into the trap. For example, using data similar to those 
in the '^''K experiments of the Zurich group, at an opti- 
cal lattice depth of 5 recoil energy for a laser wavelength 
of 825 nm, a mean trapping frequency of 80 Hz, and an 
(independently tuned) ratio \U\/t = 8, the total number 
of fermions has to be smaller than ^ 500 to have a pure 
superfiuid core. Remembering that this is only a rough 
estimate, this is nevertheless a rather small number to 



achieve under current conditions pi] . This critical num- 
ber can be larger if we allow for a central region of band 
insulator in addition to a superfiuid shell. Depending on 
the sentivity to detect the superfiuid shell, and thus the 
number of atoms that one needs to have in the superfiuid, 
the number of fermions can be larger than the above es- 
timate. The situation is better for the lighter atom ^Li. 
Using the same numbers as above, but with a laser wave- 
length of 1064 nm, the upper critical number becomes 
^ 7000, which should be feasible in current experiments. 

As for experimental signature, the observation of the 
coherence peak in the momentum distribution should sig- 
nify the onset of the BEG of pairs, and superfiuidity can 
be proven when vortices are observed when rotating the 
trap and optical lattice (although to date, vortices have 
been seen in a rotating bosonic BEG superfiuid in an opti- 
cal lattice [SS]). These are simpler probes than say, using 
noise correlation [36l [37l [38] to deduce the broken trans- 
lation symmetry of the AF state in the U > Q case|39j. 



VI. EFFECT OF DOPING 

So far, we have shown that the spin-balanced popu- 
lation for J7 < already presents several advantages to 
tackle the Mott and AF physics. But one of its main ad- 
vantages is the possibility to effectively "dope" the [/ > 
system by looking at spin imbalanced U < Q systems. 
This then allows to settle experimentally the still con- 
troversial issue of the presence or absence of the d-wave 
superfiuid (dSF) in the repulsive Hubbard model doped 
away from half-filling. As mentioned in the Introduction, 
this doping may be difficult to do directly in the U > Q 
system due to the presence of the overall harmonic trap. 
Via the transformation, the U > Q model away from half- 
filling maps to the U <Q model with an effective Zeeman 
field. In the context of cold atom experiments, this corre- 
sponds to a (fixed) imbalance of spin- up versus spin-down 
fermions, which can readily be achieved to an accuracy 
of a few percent currently |40] (see |23] for a review). 

We now examine some of the observables in that case, 
and in particular what would be the consequences of the 
existence of a d-wave superfiuid phase in the repulsive 
Hubbard model for the C/ < phase. We perform this 
analysis for the homogeneous system and will discuss the 
possible effects of the trap at the end of this section. 



A. Transformation of the operators 

Under the transformation, Eq. ([6|, the superfiuid or- 
der parameter A^"-*^ (the label a — s^d indicates the s- 
wave or d-wave symmetry of the order parameter) maps 
to the commensurate antiferromagnetic order parameter 
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stag ' 



k 

k 

k 



4^^ 



= 1 is the s-wave form factor and 
cosfcy) is the dx^-,,^ form factor. 



where previously, 

~~ (cos kx COa '^y J uxxv^ 

Thus, ij there is a regime of dSC in the U > Hubbard 
model, then correspondingly, there is a regime of dAF in 
the U < model. This is the anafogue of the well-known 
mapping at precisely half-filling of the Hubbard model 
between the ground states of sSF or sCDW at J7 < and 
the sAF at J7 > (see Fig[l]and Table |l]). 



B. Momentum distribution and noise correlation: 
a mean field calculation 

In this section, we briefly outline a simple calculation 
to illustrate experimental signatures of dAF states that 
may exist in the U < system. We need to emphasize 
from the outset that since there are no microscopic ana- 
lytical calculation of the dSF nor the dAF states in the 
2D Hubbard model, we will instead use a toy mean field 
model that does give rise to such states, in order to cal- 
culate some expected responses such as the momentum 
distribution and noise correlation [551 1571 [55] . While such 
a mean field model misses out correlations and quan- 
tum fluctuations, the objective here is to demonstrate 
that the symmetry of the CDW or SF order parameters 
has very definite signatures in noise correlation exper- 
imente[39j. Thus, we simply assume that the U < Q 
Hubbard model acquires a mean field (MF) form, which 
in momentum space becomes: 



H 



(a) 
MF 



E 



particular band structure etc. However, the following 
mean field theory only makes sense if Q = (tt/o, 7r/a) is a 
commensurate wavevector: the interaction above couples 
k with k + Q, and this in turn couples to fc -I- 2Q which 
is the same as k only for commensurate Q. 

This MF Hamiltonian is diagonalized by a Bogoliubov 
rotation: 



„t 



cos(6'fc) 4^ -I- sin(6'fe) cj'^g ^ 



Then, provided 



(13) 



cos(6'fc) 
sin(0fe) = 
^k = 



1 - 



2 flk 



idQ Mq 4>kY 



2 rife 



-,1/2 



-,1/2 



2h 



the MF Hamiltonian becomes 



Hmf = {^M^k + Eipih 



1/2 

(14) 



(15) 



with the "magnetic band" energies: 



El 



The ground state is then made up by filling these bands 
up to the respective Fermi surface /c^^. with the occupa- 
tion number = {AFQ\a\ak\AFQ) — Q{k^ — k) (and 
similarly for the /3 band): 



(16) 



crk 



(12) 



where we have introduced (by hand) the mean field order 
parameter Eq. (12 1, and as before, a labels the s or 
dx^-y^ order parameter. As usual, via a global gauge 
transformation, the order parameter can be chosen to be 
real: Mg"-* = Mq'^^ = Mq. Our main interest is in 
the phases of the U > Hubbard mode without Zeeman 
field away from half-filling. This then corresponds to the 
[/ < mean field model above (Eq. 12 1 at /i = and 
finite h. In principle, the nesting wavevector Q should 



Note that this wavefunction does not have definite num- 
ber and Ni for each species: this is a consequence of 
the quasi-particles a and (3 carrying indefinite spin, an 
analogue of the textbook number non-conserving Bogoli- 
ubov quasi-particles or the BCS wavefunction or the BCS 
mean field Hamiltonian. In turn, here for the AF, we 
have indefinite spin (but definite charge) quasi-particles 
because we have assumed the mean field decoupling to 
be in the S'^.-axis in spin space. It is straightforward to 
show (at least when /i = 0) that the same results can be 
gotten for a spin-conserving mean field wavefunction (eg. 
when the mean field decoupling is in the z spin axis). 

On the square lattice, since e^a + ^k+Qa — —2^a-, tak- 
ing the same chemical potential and bare dispersion for 



the two spin species, E'^ 



a, 13 _ 



-fj, ± fife. Hence for the 



s-wave case there is always a band-gap of size > 2gQMQ, 
and the minimum gap occurs at e/c| — Sk+Qi = 2/i. Thus 



be a variational parameter to be determined from the the half-filled lattice (1 fermion per site) with 7V| 
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has the /?-band completely filled up and an empty a- 
band: this is a magnetic insulator since adding one more 



t fermion has cj.^ ^ cos(^^fe) creating an a-particle, 
but this costs at least the band-gap energy (and similarly 
for adding a | spin particle). For the d^i _yi-viSMQ case, 
since the form factor 0^ — cos kr^ — cos fcy has nodes at 
kx = iky, the band gap also vanishes at these positions, 
giving rise to a pseudo-gap only. At less (more) than 
half-filling, the "magnetic Fermi surface" is within the 
f3- (a-) band and there is no energy gap to adding an 
extra fermion: the system is "metallic" . We should be 
careful about the meaning of E'j^'^: this is not the phys- 
ical excitation energy (unlike in the BCS MF theory): in 
particular, the spin-wave (Goldstone mode) spectrum is 
missing. 

Using the ground state Eq. ( [T6| , the T = momentum 
distribution for each spin component is: 



(4^CfcT) - cos^(0fe)< + 



sin^(6'fe) n? 



(17) 



k+Ql'^'k+Ql) 



\ek)nl+cos\ek)nl (18) 



In experiments, assuming we can image perpendicular to 
the plane of the lattice, this result has to be convolved 
with the Wannier function wik) (in momentum space) for 
the optical lattice: {nl^^^/^) \'w{k)\^ {c\^Cka) ■ (When 
there are more than 1 plane of the lattice, we also need 
to integrate over the planes.) 

The T = noise correlation is proportional to the con- 
nected correlation function: 



G„„'{k,k') = (cl^Ckacl,^,Ck'a') - {cl^Cka){cl,^,Ck'a') 

= -&a,~a'h',k+Q sin^ 6k cos^ Ok (rif (1 - )^ 



-5a ,-cr' Sk'.k+Q {QQ Mq (j)kY 



{gq Mq (f>k) 



(19) 



(Again, for the experimentally measured noise correla- 
tion, this result has to be multiphed by \w{k)\'^\w{k')\'^ .) 
We plot in Fig. [Sjthe noise correlation for the 2D dx2_y2- 
wave AF, at ii — 0. Noise correlation should clearly dis- 
tinguish between s-wave and dx2-y2 AF, thanks to the 
appearance of nodes in the ±(7r/a, tt/o) directions. 



C. Experimental observation and complications 

We thus see that the negative U side offers the great 
advantage to access directly the doped regime, without 
having to suffer directly from the presence of the trap. 
We have presented in the previous section some pos- 
sible experimental signatures that a d-wave superfluid 
phase would give when suitably transformed to the neg- 
ative side. Of course, even if the situation is poten- 
tially improved by the transformation to the negative U 




FIG. 5: Density (grey-scale) plot of the noise correlation for 
the 2D _y2-wa,ve AF against kx and ky when k' = k + Q. 
Parameters used: qqMq = 0.10, n = 0, h = 0.23. 



side, probing the doped phase is a considerable challenge. 
Some of the limitations are obvious. The most immedi- 
ate one is the effect of the temperature. Indeed, already 
for the antiferromagnetic phase, the temperature is in 
competition with an energy of order J ~ it'^ /U the ki- 
netic energy of the pairs. When looking at the doping 
effects one has to face even smaller energies. Lowering 
the temperature is thus a must. 

The second limitation is again the trap. Although the 
mapping to the negative side avoids the direct effect of 
the trap on the doped holes, the trap still has a poten- 
tially indirect effect. Indeed as we saw in the section 
on antiferromagnetism, the trap will act, for the U > 
side (model Rl), as an effective Zeeman field, and lead to 
two shell regions, with fully polarized spins up or down 
(Fig. [4 1. When holes are introduced into the system of 
model Rl, they will thus have the possibility to go in one 
of these two polarized regions or in the antiferromagnetic 
region, which corresponds to the shell where the effective 
Zeeman field is not strong enough to fully polarize the 
system. Where the holes go will depend on how much 
kinetic and interaction energy they can gain in the three 
different phases since they are not sensitive directly to 
the presence of the Zeeman field. A naive calculation 
neglecting the presence of the energy scale J, leads di- 
rectly to the holes going in the interface between the 
fully polarized up and fully polarized down regions. In 
the U < language, i.e. for model Al, this corresponds 
to the excess of one spin species going to the edge of the 
band insulator region. This seems to suggest, that when 
the energy scale J is put back in the problem, the holes 
will indeed go into the antiferromagnetic region, and thus 
can lead potentially to the d-wave superconducting phase 
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there. This is however a delicate question since the ki- 
netic energy of a hole in a ferromagnetic environment is 
in principle higher due to the absence of frustration of 
the antiferromagnetic order upon hole motion. This im- 
portant question thus fully deserves more analytical and 
numerical study. It remains however academic until se- 
rious progress on the temperature issue have been made. 

VII. SUMMARY, CONCLUSION 

In this article, we have explored in depth the possibil- 
ities offered by quantum simulating the attractive model 
(Al) in cold atoms in optical lattices, to understand the 
physics of the repulsive Hubbard model, via a well known 
canonical particle-hole transformation. We have argued 
that there are certain advantages in doing experiments 
in the attractive regime. 

For the undoped case the attractive side replaces the 
Mott phase and the antiferromagnetic phase by a phase 
composed only of pairs (Mott insulator for U > 0) and 
that undergoes a superfluid transition (antiferromagnet 
for U > 0). The trap which exists in any realistic ex- 
periment does not really affect the observation of these 
two phases since it can only add a core of band insulator 
at the center, thereby not spoiling the observation of the 
superfluid. The attractive side also offers the advantage 
of only having to test for the pairing for the observation 
of the band insulator that must be simpler than testing 
for the absence of doubly occupied and empty sites on 
the repulsive side. 

Another key advantage of the attractive side is the 
relative ease in controlling the spin population imbal- 
ance. Via the canonical transformation, this corresponds 



to doping away from half filling for the repulsive side, 
which is hard to achieve because of the presence of the 
harmonic trapping potential that moves the holes away 
from the central region of the trap. This doping in the 
repulsive regime is needed to quantum simulate and ex- 
plore the possibility of a d-wave superfluid, the funda- 
mental question of great relevance to the cuprate high 
temperature superconductors. This question can be an- 
swered instead in the attractive regime by exploring the 
presence or absence of a d-wave antiferromagnet. We 
indicate in this paper several ways to probe for the exis- 
tence of such a phase. We have also pointed out a number 
of technical difficulties of the proposed approach and, in 
some cases, suggest possible solutions. 
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